/*! 
 ***************************************************************************
 * \file FiveSpot3D.cpp
 * Two phase flow in porous media in 3D.
 *
 * \par Description
 * This code solves the two phase flow in porous media in 3D, the
 * equations are posed in pressure-saturation formulation:
 * 
 * These equation are solved using the IMPES algorithm.
 * The domain of study is 
 * \f$ (x,y,z) \in [0,182.76 m] \times [0,182.76 m]\times[0,9.14 m] \f $, and is
 * initially filled with oil. Then, water is injected on a corner and
 * the oil is pushed to the opposed corner of the domain.
 *
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./FiveSpot3D  \endverbatim
 ***************************************************************************
 *  \author    Luis M. de la Cruz [ May 21 10:34:21 GMT 2011 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/TwoPhaseEquation.hpp"
#include "./FS/FSIP1.hpp"
#include "./FS/FSES1.hpp"
#include "Solvers/TDMA.hpp"
#include "Solvers/explicit.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in two dimensions
typedef TunaArray<double,3>::huge ScalarField3D;

timer crono;

int pausa;

int main() 
{
  double length_x, length_y, length_z, dt, Tmax, tolerance = 1e-6; 
  int num_nodes_x, num_nodes_y, num_nodes_z, frequency, tdma_iter = 20;
  double porosity, permeability, mu_w, mu_o, Srw, Sro;
  double injection, pres_out, Sw_in;
  int idt_p;  

  string filenameInput;

  // Reading data from the input file
  //  cout << "\n ---> File name for input: "; cin >> filenameInput;
  //  std::ifstream input_cin (filenameInput.c_str()); 
  std::ifstream input_cin ("inputFS3D"); 
  input_cin >> length_x >> length_y >> length_z
	    >> num_nodes_x >> num_nodes_y >> num_nodes_z
	    >> porosity
	    >> permeability
	    >> mu_w
	    >> mu_o
	    >> Srw
	    >> Sro
	    >> injection
	    >> pres_out
	    >> Sw_in
	    >> dt
	    >> Tmax
	    >> idt_p
	    >> frequency;
  input_cin.close();
  
    double Gamma = permeability;

    StructuredMesh<Uniform<double, 3> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y,
					     length_z, num_nodes_z);
    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    double dz = mesh.getDelta(Z);
    mesh.print();

    ScalarField3D p ( mesh.getExtentVolumes() );
    ScalarField3D Sw ( mesh.getExtentVolumes() ); 
    ScalarField3D p_n( mesh.getExtentNodes()   );
    ScalarField3D Sw_n( mesh.getExtentNodes()   );

    // Initial conditions
    Sw = 0;
    p = pres_out; 
 
    Range all = Range::all();

    NumUtils::interpolateToNodes(p_n, p);
    InOut::writeToFile_DX(p_n, 0, "./DataFS3D/pres.", dx, dy, dz);
    NumUtils::interpolateToNodes(Sw_n, Sw);
    InOut::writeToFile_DX(Sw_n, 0, "./DataFS3D/satu.", dx, dy, dz);
        
    SparseMatrix< Diagonal<double, 3> > A(num_nodes_x, num_nodes_y, num_nodes_z); 
    ScalarField3D                       b(num_nodes_x, num_nodes_y, num_nodes_z);

    TwoPhaseEquation< FSIP1<double, 3> > pressure(p, A, b, mesh.getDeltas());
    pressure.setDeltaTime(dt); 
    pressure.setPermeability(permeability);
    pressure.setPorosity(porosity);
    pressure.setSrw(Srw);
    pressure.setSro(Sro);
    pressure.setViscosity_w(mu_w);
    pressure.setViscosity_o(mu_o);
    pressure.setInjection(injection);
    pressure.setNeumann  (LEFT_WALL);
    pressure.setNeumann  (RIGHT_WALL);
    pressure.setNeumann  (TOP_WALL);
    pressure.setNeumann  (BOTTOM_WALL);
    pressure.setNeumann  (FRONT_WALL);
    pressure.setNeumann  (BACK_WALL);
    pressure.setSaturation(Sw);
    pressure.print();

    //
    // In these case, the boundary conditions are calculated in
    // the numerical scheme (the adaptor BLES7).
    //
    TwoPhaseEquation< FSES1<double, 3> > saturation(Sw, A, b, mesh.getDeltas());
    saturation.setDeltaTime(dt); 
    saturation.setPermeability(permeability);
    saturation.setPorosity(porosity);
    saturation.setSrw(Srw);
    saturation.setSro(Sro);
    saturation.setInjection(injection);
    saturation.setViscosity_w(mu_w);  
    saturation.setViscosity_o(mu_o); 
    saturation.setNeumann  (LEFT_WALL);
    saturation.setNeumann  (RIGHT_WALL);
    saturation.setNeumann  (TOP_WALL);
    saturation.setNeumann  (BOTTOM_WALL);    
    saturation.setNeumann  (FRONT_WALL);
    saturation.setNeumann  (BACK_WALL);
    saturation.setPressure(p);
    saturation.print();


    ofstream error_fp("./DataFS3D/error"), 
      residual_fp("./DataFS3D/residual");
    int iteration = 1, idprint = 1;
    double t = dt;
    double error = 1.0;
    double residual = 1.0;

    int pausa;
    double maximo;

    int bi = saturation.get_bi(), bj = saturation.get_bj(), bk = saturation.get_bk();
    int ei = saturation.get_ei(), ej = saturation.get_ej(), ek = saturation.get_ek(); 

    frequency *= (86400 / dt);

    double total_time = 0, pres_time = 0, sat_time = 0;
    while (t <= Tmax) {

      if ( !(iteration % idt_p) ) {
	crono.tic();
	pressure.calcCoefficients();
	Solver::TDMA3D(pressure, tolerance, tdma_iter, 1.0);
	pres_time += crono.toc();
	
	residual = pressure.calcResidual();
	error = pressure.calcErrorL1();
	
	error_fp << iteration << "\t" << error << std::endl;
	residual_fp << iteration << "\t" << residual << std::endl;	
	
	pressure.update();
      }  

      crono.tic();
      saturation.calcCoefficients();
      Solver::solExplicit3D(saturation);
      sat_time += crono.toc();

      error = saturation.calcErrorL1();
      saturation.update();


// ----- Flux-limiter to cut it off the wiggles
      /*
      // Injector treatment: Extrapolation
      Sw(bi,bj,bk) =  2 * Sw(bi+1, bj+1, bk) - Sw(bi+2, bj+2, bk);
      Sw(bi,bj,bk+1) = Sw(ei, ej, ek);
      Sw(bi,bj,bk-1) = Sw(ei, ej, ek);

      Sw(ei+1, ej+1, ek) = Sw(ei, ej, ek);
      Sw(ei+1,ej+1,ek+1) = Sw(ei, ej, ek);
      Sw(ei+1,ej+1,ek-1) = Sw(ei, ej, ek);
      /* */
      // Producer treatment: Extrapolation
      Sw(ei,ej,ek) =  2 * Sw(ei-1, ej-1, ek) - Sw(ei-2, ej-2, ek);
      Sw(ei,ej,ek+1) = Sw(ei, ej, ek);
      Sw(ei,ej,ek-1) = Sw(ei, ej, ek);

      Sw(ei+1, ej+1, ek) = Sw(ei, ej, ek);
      Sw(ei+1,ej+1,ek+1) = Sw(ei, ej, ek);
      Sw(ei+1,ej+1,ek-1) = Sw(ei, ej, ek);

      for(int kk = 0; kk <= (ek+1); ++kk)
	for(int j = 0; j <= (ej+1); ++j)
	  for(int i = 0; i <= (ei+1); ++i) {
	    if ( Sw(i,j,kk) > Sw_in) Sw(i,j,kk) = Sw_in;
	    if ( Sw(i,j,kk) < Srw ) Sw(i,j,kk) = Srw;
	  }
// ----- 

      saturation.updatePhi(Sw);

      
      if ( !(iteration % frequency) ) {
	/*
	NumUtils::interpolateToNodes(p_n, p);
	InOut::writeToFile_DX(p_n, idprint, "./DataFS3D/pres.",dx, dy, dz);
	
	NumUtils::interpolateToNodes(Sw_n, Sw);
	InOut::writeToFile_DX(Sw_n, idprint, "./DataFS3D/satu.",dx, dy, dz);
	*/

	InOut::writeToFile_DX(p, idprint, "./DataFS3D/pres.",dx, dy, dz);
	InOut::writeToFile_DX(Sw, idprint, "./DataFS3D/satu.",dx, dy, dz);

	std::cout << "\n\n ---- ID print => " << idprint
		  << "\t - Iteration => " << iteration
		  << "\t Step = " << t
		  << "\n ---- Res = " << residual
		  << "\t Err = " << error;
	std::cout.flush();
	++idprint;
      }

      ++iteration;
      t += dt;
    }

    cout << "\n\n Presure Elapsed time = " << pres_time << " segs \n";
    cout << "\n Saturation Elapsed time = " << sat_time << " segs \n";
    
    std::cout << "\n\n NOTE: The data were written in the [./DataFS3D] directory \n\n";

    return 0;
}




